###
# Replication for Simmons, Lloyd and Stewart
###


# A Note on the Data: 
# The data was collected to run through criminalizations in 2008
# we use a particular hack of listing the criminalization date
# of countries that we know have not yet criminalized as 2012 and
# list as NA those whose criminalization status is unknown entirely.
# This works for this code but should obviously be addressed if
# one is attempting to update the data.

library(survival)
library(foreign)
library(data.table)
library(Greg)

# A function to make tables we can use in appendix or for the paper
# simply calculates hazard ratios of a cox model along with CI's and 
# p-values.  Either returns the table or the xtable 
returnTable <- function(mod,rownames=NULL,xtable=FALSE, digits=3) {
  tab <- summary(mod)$coefficients[,c(1,2,6)]
  Lower <- exp(tab[,1] + qnorm(.025)*tab[,2])
  Upper <- exp(tab[,1] + qnorm(.975)*tab[,2])
  HR <- exp(tab[,1])
  tab <- data.frame(HR, Lower, Upper, pvalue=tab[,3])
  rownames(tab) <- rownames
  if(xtable) {
    return(xtable::xtable(tab, digits=digits))
  } else {
    return(round(tab, digits))
  }
}

#prints out number of observations, number of events,
#and the cox non-proportional hazards tests of Grambsch and Therneau (1994)
#for a cox model
diagnose <- function(mod) {
  n <- summary(mod)$n
  ne <- summary(mod)$nevent
  p.nph <- cox.zph(mod)$table[,3] #grab p-values
  global <- p.nph[length(p.nph)] < .05
  nph <- any(p.nph < .05)
  cat(sprintf("N Complete Cases: \t %i \nN Events: \t \t %i \nAny NPH Prop p<.05: \t %s \nGlobal NPH Prop p<.05: \t %s",
              n, ne, nph, global))
}

#Make a little function that will add p-values
#to an ordinal probit table created by MASS.  Return
#a rounded table out to 3 digits.  Also provides number of obs
op_summary <- function(mod) {
  m.coef <- data.frame(coef(summary(mod)))
  m.coef$pval <- round((pnorm(abs(m.coef$t.value), lower.tail = FALSE) * 2),2)
  list(coef=round(m.coef,3), n=mod$n)
}

make_na_data <- function(form, data) {
  mini <- data[,all.vars(as.formula(form)),with=FALSE]
  na.omit(mini)
}


###
# Figure 1
###
# here we use updated data reflecting retrospective 
# criminalization information. Thus it differs from
# the data available at the time.
data <- data.table::fread("updatedcrimdata.csv")
colnames(data) <- c("name", "crimdate")
x <- cumsum(table(data[, mean(crimdate),by=name]$V1))
x <- x[which(names(x)==1997):length(x)]
library(ggplot2)
x <- data.frame(year=1997:2015, y=x)
pdf("../figures/Figure1.pdf", width=6, height=5)
ggplot(x, aes(x=year,y=y)) + geom_line() + geom_point() + theme_bw() +
  ylab("Cumulative Countries Criminalized") + xlab("Year") + 
  scale_x_continuous(breaks=seq(1997,2015,by=1)) + 
  scale_y_continuous(breaks=seq(0, 145, by=10)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
dev.off()



########
# Data Cleaning
########

load("CompleteDataWithWeights.RData")
data <- data.table(data)


##
#Create the survival analysis variables
##

#Start timer at 1990
data$t <- data$year - 1990
data$t0 <- data$year - 1991
subset <- data[year>=1990,]
setkey(subset, country, year)

#generate lags and fail conditions
subset[,crim_ht_lag1 := shift(crim_ht), by=country]
subset[,fail := ifelse(crim_ht_lag1==0 & crim_ht==1, 1, 0)]
subset[,drop := ifelse(crim_ht_lag1==1, 1, 0)]
#drop out those who have already failed
subset <- subset[drop!=1,]
subset$drop <- NULL
#rename
data <- subset

save(data, file="HTreadytomodel.RData")

#Figure 2 can be replicated by the URL:
# https://books.google.com/ngrams/graph?content=transnational+crime%2Ctransnational+organized+crime&year_start=1970&year_end=2008&corpus=15&smoothing=3&share=&direct_url=t1%3B%2Ctransnational%20crime%3B%2Cc0%3B.t1%3B%2Ctransnational%20organized%20crime%3B%2Cc0

devtools::install_github("ngramr", "seancarmody")
require(ggplot2)
pdf("../figures/Figure2.pdf", width=8, height=4)
ng  <- ngramr::ngram(c("transnational crime","transnational organized crime"), year_start = 1970,year_end=2010)
ggplot(ng, aes(x=Year, y=Frequency, colour=Phrase, linetype=Phrase)) +
  geom_line() + theme_bw() + scale_x_continuous(breaks=seq(1970,2010,by=10))
dev.off()


########
# Basic models
########

#load("HTreadytomodel.RData")

#########
##Table 1, Model 1
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2) + uspressure_1 + 
         pspline(rule_of_law, df=3) + rattocprot + cluster(country)
t1m1 <- coxph(form, data=data, x=TRUE)
tab <- returnTable(t1m1, xtable=FALSE,
                   rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                              "US Pressure", "Rule of Law",
                              "Rule of Law Nonlinear", "Ratification"))
tab
diagnose(t1m1)

##
# Figure 3
##
pdf(file="../figures/Figure3.pdf", width=10, height=4)
par(mfrow=c(1,2))
Greg::plotHR(t1m1, xlim=c(0,60), ylim=c(.5, 6), rug="ticks", 
             main="Spatial Spline", 
             xlab="Number of Roads to Criminalized Neighbors",
             col.term="black",col.se="lightgrey")
Greg::plotHR(t1m1, xlim=c(0,10), ylim=c(1, 2.5), rug="ticks", 
             main="Spatial Spline (Zoomed)", 
             xlab="Number of Roads to Criminalized Neighbors",
             col.term="black",col.se="lightgrey")
dev.off()

#########
#Table 1, Model 2
form <- Surv(t0,t,fail) ~ pspline(Crimsum, df=2) + uspressure_1 + pspline(rule_of_law, df=3) + 
                          rattocprot + pspline(USaid_gdp, df=2) +
                          useimfcr + 
                          pspline(us_trade_share_cap, df=2) + 
                          pspline(eu_trade_share_cap, df=2) +
                          cluster(country)
t1m2 <- coxph(form, data=data, x=TRUE)

tab <- returnTable(t1m2, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear", "Ratification",
                                    "US aid GDP", "US Aid Nonlinear",
                                    "IMF Credit", "US Trade Share", 
                                    "US Trade Share Nonlinear",
                                    "EU Trade Share", "EU Trade Share Nonlinear"),
                   xtable=FALSE)
tab

diagnose(t1m2)

#######
#Table 1, Model 3

form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2) + uspressure_1 + 
                          pspline(rule_of_law, df=3) + rattocprot + 
                          pspline(ichild_labor, df=3) + mid + 
                          pspline(remittances,df=3) +
                          cluster(country)

t1m3 <- coxph(form, data=data, x=TRUE)

tab <- returnTable(t1m3, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                             "US Pressure", "Rule of Law",
                             "Rule of Law Nonlinear", "Ratification",
                             "Child Labor", "Child Labor Nonlinear",
                             "Middle Income Country",
                             "Remittances", "Remittances Nonlinear"),
                   xtable=FALSE)
tab

diagnose(t1m3)

#########
#Table 1, Model 4
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2) + uspressure_1 + 
                          pspline(rule_of_law, df=3) + rattocprot + 
                          pspline(ichild_labor,df=3) + mid +
                          pspline(par_women_lower_, df=3) +
                          cluster(country)
t1m4 <- coxph(form, data=data, x=TRUE)

tab <- returnTable(t1m4, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                             "US Pressure", "Rule of Law",
                             "Rule of Law Nonlinear", "Ratification",
                             "Child Labor", "Child Labor Nonlinear",
                             "Middle Income Country",
                             "Share Women in Parliament",
                             "Share Women Nonlinear"),
                   xtable=FALSE)
tab
diagnose(t1m4)




#########
#Table 1, Model 5
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2) + uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot+ 
  pspline(latentmean, df=3) +
  cluster(country)

t1m5 <- coxph(form, data=data, x=TRUE)

tab <- returnTable(t1m5, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                             "US Pressure", "Rule of Law",
                             "Rule of Law Nonlinear", "Ratification",
                             "Respect For Human Rights",
                             "Respect for HR Nonlinear"),
                   xtable=FALSE)
tab

diagnose(t1m5)

############
#Table 1, Model 6
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2) + uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + 
  pspline(crimsumany_HRIO, df=3) +
  cluster(country)
t1m6 <- coxph(form, data=data, x=TRUE)

tab <- returnTable(t1m6, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                             "US Pressure", "Rule of Law",
                             "Rule of Law Nonlinear", "Ratification",
                             "Shared HROrg Crim",
                             "Shared HROrg Crim Nonlinear"),
                   xtable=FALSE)
tab

diagnose(t1m6) #global test is fine, an individual test rejects
#where is the failure happening?
cox.zph(t1m6) #it is parts of the spline for the shared HR orgs. 

##########
# Table 2: Robustness of Diffusion via Roads: Other "Proximity" Measures
##########

#Model 1: Policy Weighted by Trade Partners
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + pspline(crim_ht_w1, df=2) +
  cluster(country)
t2m1 <- coxph(form, data=data, x=TRUE)

tab <- returnTable(t2m1, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear", 
                                    "Policy Weighted by Trade Partners",
                                    "Policy Weighted by Trade Nonlinear"),
                   xtable=FALSE)
tab

diagnose(t2m1)

#Model 2: Regional Press Stories on HT
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + pspline(laglogreg_arts2, df=2) +
  cluster(country)
t2m2 <- coxph(form, data=data, x=TRUE)

tab <- returnTable(t2m2, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear", 
                                    "Regional News",
                                    "Regional News Nonlinear"),
                   xtable=FALSE)
tab

diagnose(t2m2)

#Model 3: Developmental Emulation
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + pspline(devlevDens, df=2) + 
  cluster(country)
t2m3 <- coxph(form, data=data, x=TRUE)


tab <- returnTable(t2m3, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear", 
                                    "Density in DevLevel",
                                    "Density in DevLevel Nonlinear"),
                   xtable=FALSE)
tab

diagnose(t2m3)

#Model 4: Civilizational Emulation
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + pspline(civDens, df=2)  + cluster(country)
t2m4 <- coxph(form, data=data, x=TRUE)

tab <- returnTable(t2m4, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                               "US Pressure", "Rule of Law",
                             "Rule of Law Nonlinear",
                             "Civilizational Crim",
                             "CivCrim Nonlinear"),
                   xtable=FALSE)
tab
diagnose(t2m4)

####
#Model 5: Legal Families
####
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + pspline(legalDens, df=2)  + cluster(country)
t2m5 <- coxph(form, data=data, x=TRUE)

tab <- returnTable(t2m5, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                             "US Pressure", "Rule of Law",
                             "Rule of Law Nonlinear",
                             "Legal Crim",
                             "LegCrim Nonlinear"),
                   xtable=FALSE)
tab

diagnose(t2m5)

##########
# Table 3 (Model 1) / Figure 4: Robustness to Other Outcomes
##########

# DV: Prosecution (ordered probit)
form <- factor(prosecution) ~ asinh(Crimsum) + rule_of_law + hr_treaties + 
        par_women_lower_ + laglogcrime_arts2 + laglogvictim_arts2
mod1 <- MASS::polr(form, data=data, method="probit")
op_summary(mod1)

# DV: Protection (ordered probit)
form <- factor(protection) ~ asinh(Crimsum) + rule_of_law + hr_treaties + 
  par_women_lower_ + laglogcrime_arts2 + laglogvictim_arts2
mod2 <- MASS::polr(form, data=data, method="probit")
op_summary(mod2)

#Rerun the models with zelig in order to make the plots
#we have to basically trick Zelig into doing the transformation
library(zeligverse)
set.seed(08540)
data$Spatial <- asinh(data$Crimsum)
xlowvalue <- 0 #set range of values to consider
xhighvalue <- 5
form <- factor(prosecution) ~ Spatial + rule_of_law + hr_treaties + 
  par_women_lower_ + laglogcrime_arts2 + laglogvictim_arts2
z.out <- zelig(form, model="oprobit", data=data)
x.low <- setx(z.out, Spatial=asinh(xlowvalue)) #asinh(0) evaluates to 0, just for clarity
x.high <- setx(z.out, Spatial=asinh(xhighvalue))
#do a lot of simulations as we are going to calculate by quantiles
s.out <- sim(z.out, x = x.low, x1 = x.high,num = 25000)
#Zelig doesn't seem to do risk ratios any more
#we want to calculate probability of high / probability of low
rr <- s.out$sim.out$x1$ev[[1]]/s.out$sim.out$x$ev[[1]]
#store values for plotting
pros_calc <- rbind(colMeans(rr),apply(rr, 2, quantile, c(.025, .975)))

#now repeat for protection
form <- factor(protection) ~ Spatial + rule_of_law + hr_treaties + 
  par_women_lower_ + laglogcrime_arts2 + laglogvictim_arts2
z.out <- zelig(form, model="oprobit", data=data)
x.low <- setx(z.out, Spatial=asinh(xlowvalue))
x.high <- setx(z.out, Spatial=asinh(xhighvalue))
#do a lot of simulations as we are going to calculate by quantiles
s.out <- sim(z.out, x = x.low, x1 = x.high,num = 25000)
#Zelig doesn't seem to do risk ratios any more
#we want to calculate probability of high / probability of low
rr <- s.out$sim.out$x1$ev[[1]]/s.out$sim.out$x$ev[[1]]
#store values for plotting
prot_calc <- rbind(colMeans(rr),apply(rr, 2, quantile, c(.025, .975)))

# get rid of that spatial variable we created
data$Spatial <- NULL
##
# Make the actual Plot
##
pdf(file = "../figures/Figure4.pdf", width = 6, height = 5) 
par(mfrow=c(1,1))

plot(seq(1, by=3, length=5), log(pros_calc[1,]), xlab="Level of Prosecution/Protection", ylab="Risk Ratio for 0 to 5 road increase", xaxt="n", yaxt="n", xlim=c(0, 14), ylim=c(-2, log(6)), main="Diffusion for Prosecution/Protection")
segments(seq(1, by=3, length=5), log(pros_calc[2,]),seq(1, by=3, length=5),log(pros_calc[3,]), lty=1)
axis(1, at=seq(1.5, by=3, length=5),labels=c("Weakest (1)", "2", "3", "4", "Strongest (5)"))
axis(side = 2 , at = axTicks(2) , label = round( exp(axTicks(2)) , digits = 3))

points(seq(2, by=3, length=5), log(prot_calc[1,])) 
segments(seq(2, by=3, length=5), log(prot_calc[2,]),seq(2, by=3, length=5),log(prot_calc[3,]), lty=2)
abline(h=0, lty=2)
legend("topleft", c("Prosecution (95% CI)", "Protection (95% CI)"),lty=c(1,2)) 
dev.off()

#########
# Table 3 (Model 2) / Figure 5: HT and ML Comparison
#########
##
#Human Trafficking DV
##
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ 
  pspline(rule_of_law, df=3) + devlev  + cluster(country)
t3m2a <- coxph(form, data=data, x=TRUE)

tab <- returnTable(t3m2a, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "Rule of Law",
                                    "Rule of Law Nonlinear",
                                    "DevLev"),
                   xtable=FALSE)
tab

diagnose(t3m2a)

##
#Money laundering 
##
# we now need to load a new dataset and recreate
# the survival model variables

htdata <- data
load("CompleteDataWithWeights.RData")
data <- data.table(data)

#Start timer at 1990
data$t <- data$year - 1990
data$t0 <- data$year - 1991
subset <- data[year>=1990,]
setkey(subset, country, year)

#generate lags and fail conditions
subset[,crim_ml_lag1 := shift(ml_crim_lax), by=country]
subset[,fail := ifelse(crim_ml_lag1==0 & ml_crim_lax==1, 1, 0)]
subset[,drop := ifelse(crim_ml_lag1==1, 1, 0)]
#drop out those who have already failed
subset <- subset[drop!=1,]
subset$drop <- NULL
#rename
MLdata <- subset
data <- htdata

form <- Surv(t0,t,fail) ~ pspline(MLsum,df=2)+ 
  pspline(rule_of_law, df=3)+ devlev + cluster(country)
t3m2b <- coxph(form, data=MLdata, x=TRUE)
tab <- returnTable(t3m2b, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                     "Rule of Law",
                                     "Rule of Law Nonlinear",
                                     "DevLev"),
                   xtable=FALSE)
tab

diagnose(t3m2b)


##
# Figure 5
##
pdf(file="../figures/Figure5.pdf", width=10, height=4)
par(mfrow=c(1,2))
Greg::plotHR(t3m2b, xlim=c(0,60), ylim=c(.5, 12), rug="ticks", 
             main="Money Laundering", 
             xlab="Number of Roads to Criminalized Neighbors",
             col.term="black",col.se="lightgrey")

Greg::plotHR(t3m2a, xlim=c(0,60), ylim=c(.5,12), rug="ticks", 
             main="Human Trafficking", 
             xlab="Number of Roads to Criminalized Neighbors",
             col.term="black",col.se="lightgrey")
dev.off()


########
# Table 4/Figure 6: By Exposure Category
########

#now we are back to human trafficking so we load back up
# that data
load("HTreadytomodel.RData")

##
#Destination Countries
##
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + cluster(country)
t4m1 <- coxph(form, data=data[data$destination==1,], x=TRUE)

tab <- returnTable(t4m1, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear",
                                    "Ratification"),
                   xtable=FALSE)
tab

diagnose(t4m1)

##
#Origin Countries
##
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + cluster(country)
t4m2 <- coxph(form, data=data[data$origin==1,], x=TRUE)

tab <- returnTable(t4m2, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear",
                                    "Ratification"),
                   xtable=FALSE)
tab
diagnose(t4m2)


##
#Transit Countries
##
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + cluster(country)
t4m3 <- coxph(form, data=data[data$transit==1,], x=TRUE)

tab <- returnTable(t4m3, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear",
                                    "Ratification"),
                   xtable=FALSE)
tab

diagnose(t4m3)


##
#internal
##
form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2)+ uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + cluster(country)
t4m4 <- coxph(form, data=data[data$internal_trafficking==1,], x=TRUE)

tab <- returnTable(t4m4, rownames=c("Sum of Roads", "Sum of Roads Nonlinear",
                                    "US Pressure", "Rule of Law",
                                    "Rule of Law Nonlinear",
                                    "Ratification"),
                   xtable=FALSE)
tab

diagnose(t4m4)


pdf(file="../figures/Figure6.pdf", width=10, height=4)
par(mfrow=c(1,2))
Greg::plotHR(t4m1, xlim=c(0,60), ylim=c(.5, 12), rug="ticks", 
             main="Destination", 
             xlab="Number of Roads to Criminalized Neighbors",
             col.term="black",col.se="lightgrey")
# 
# Greg::plotHR(t4m2, xlim=c(0,60), ylim=c(.5, 12), rug="ticks", 
#              main="Origin", 
#              xlab="Number of Roads to Criminalized Neighbors",
#              col.term="black",col.se="lightgrey")

Greg::plotHR(t4m3, xlim=c(0,60), ylim=c(.5, 12), rug="ticks", 
             main="Transit", 
             xlab="Number of Roads to Criminalized Neighbors",
             col.term="black",col.se="lightgrey")
# 
# Greg::plotHR(t4m4, xlim=c(0,60), ylim=c(.5, 12), rug="ticks", 
#              main="Internal Trafficking", 
#              xlab="Number of Roads to Criminalized Neighbors",
#              col.term="black",col.se="lightgrey")
dev.off()


########
# Figure 7: Weighting neighbors by type.
########

modeld <- coxph(Surv(t0,t,fail) ~ pspline(originSum, df=2) + 
                  pspline(destinationSum, df=2)+ 
                  transitSum +
                  pspline(rule_of_law, df=3) +
                  uspressure_1 + rattocprot  + 
                  cluster(country), data=data)

tab <- returnTable(modeld, rownames=c("Origin", "Origin Nonlinear",
                                    "Destination", "Destination Nonlinear",
                                    "Transit", 
                                    "Rule of Law", "Rule of Law Nonlinear",
                                    "US Pressure", "Ratification"),
                   xtable=FALSE)
tab

pdf(file = "../figures/Figure7.pdf", width = 10, height = 4) 
par(mfrow=c(1,2))
Greg::plotHR(modeld, xlim=c(0,60), ylim=c(.5, 8), rug="ticks", 
             main="Origin Neighbors", 
             xlab="Number of Roads to Criminalized Origin Neighbors",
             col.term="black",col.se="lightgrey")
Greg::plotHR(modeld, term=2, xlim=c(0,60), ylim=c(.5, 8), rug="ticks",
             main="Destination Neighbors",
             xlab="Number of Roads to Criminalized Destination Neighbors",
             col.term="black",col.se="lightgrey")
dev.off()


########
# Claims in the paper
########

#Model 4: country with no women in parliament is about 10 times less likely
# to criminalize than median country 11.1 women.
# also compare median to amost 50.


form <- Surv(t0,t,fail) ~ pspline(Crimsum,df=2) + uspressure_1 + 
  pspline(rule_of_law, df=3) + rattocprot + 
  pspline(ichild_labor,df=3) + mid +
  pspline(par_women_lower_, df=3) +
  cluster(country)
mini <- make_na_data(form, data)
mod <- coxph(form, data=mini, x=TRUE)
temp <- data.frame(rbind(mini[1,],mini[1,]))
temp$par_women_lower_ <- c(10,48)
pred <- predict(mod, temp, type="risk")
pred[2]/pred[1]

# Get Rid of the temporary file we created
unlink("HTreadytomodel.RData")

# Some code to help in plotting hazard ratios when there is missingness.
# mini <- make_na_data(form, data)
# mod <- coxph(form, data=mini, x=TRUE)
# Greg::plotHR(mod, term=4, xlim=c(0, 1), rug="ticks",
#              main="DV",
#              xlab="Number of Roads to Criminalized Destination Neighbors",
#              col.term="black",col.se="lightgrey")

